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Abstract. Reaction diffusion systems are often used to study pattern formation in 
biological systems. However, most methods for understanding their behavior are chal- 
lenging and can rarely be applied to complex systems common in biological applications. 
I present a relatively simple and efficient, non-linear stability technique that greatly aids 
such analysis when rates of diffusion are substantially different. This technique reduces a 
system of reaction diffusion equations to a system of ordinary differential equations track- 
ing the evolution of a large amplitude, spatially localized perturbation of a homogeneous 
steady state. Stability properties of this system, determined using standard bifurcation 
techniques and software, describe both linear and non-linear patterning regimes of the 
reaction diffusion system. I describe the class of systems this method can be applied 
to and demonstrate its application. Analysis of Schnakcnborg and substrate inhibition 
models is performed to demonstrate the methods capabilities in simplified settings and 
show that even these simple models have non-linear patterning regimes not previously 
detected. Analysis of a protein regulatory network related to chemotaxis shows its ap- 
plication in a more complex setting where other non-linear methods become intractable. 
Predictions of this method are verified against results of numerical simulation, linear 
stability, and full PDE bifurcation analyses. 



1. Introduction 

Reaction diffusion equations (RDEs) have provided a ubiquitous framework for studying 
pattern formation in chemical and biological systems [40l [101 EH 1321 EQI [H] . As a result of 
their sustained interest, numerous linear [40], weakly non-linear [SH [39l [36l [21] , and fully 
non-linear [171 [Ml [Ml [H] [71 [23l [30] techniques for analyzing RDEs have been developed. 
Here I present an efficient, relatively simple addition to this toolbox, aimed at analyzing 
systems where diffusivities are substantially different. 

A difference in rates of diffusion has been implicated as being of vital importance for 
patterning in numerous biological systems. In the context of cell biology, some rates of 
diffusion are not just different but are vastly different, varying by factors of 100 — 1000. 
Numerous cell functions arc controlled by regulators, such as "GTPase's" in cell motility 
[20l[Tl[4l[28l[13l[29], "ROP's" in plant development [9], and "Min" proteins in bacterial 
division [151 116j . All of these regulatory proteins have fast and slow diffusing components 
since they exist in membrane bound and unbound states. Analysis of these types of systems 
motivated the development of this method. 
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Here I consider a generic system of RDEs with a large diffusion disparity and highlight a 
useful method for understanding their linear and non-linear stability properties. Consider 

du 

(1) -Q^{x,t) ^ f{u,v;p) + T>^u^^, 

where u and v are vectors of slow and fast diffusing variables respectively, Du,Dv are 
diagonal matrices of diffusion coefficients, and p is a vector of reaction parameters. 

The "Local Perturbation Analysis" (LPA), is a non-linear stability technique applicable 
to systems of this type. This method, originally devised by Marcc and Grieneisen [T2], is 
a bridge between linear and non-linear analysis methods having benefits of each. Linear 
stability analysis (LSA) j40| is straitforward and widely used, but is limited to providing 
linear information. Non-linear methods, while more informative, are much more challenging, 
specific to the system being investigated, and rarely scale up to complex systems with many 
variables. The LPA provides non-linear stability information beyond that of LSA, but is 
relatively simple to implement. 

In contrast to LSA which probes stability of a homogeneous steady state (HSS) with 
respect a small amplitude, spatially extended perturbation, the LPA probes stability with 
respect to a spatially localized, large amplitude perturbation of the slow variable u. As 
will be shown, when diffusion of u (resp. v) is sufRcicntly slow (rcsp. fast) the perturbed 
region and broader domain evolve according to an approximate collection of ODE's on the 
timescale of reactions 

(2a) ^i^^^t)^ f{uay-p), 

(2b) —^^^t)=giuS,vS;p), 
(2c) ^^^^t) = f{u\va;p). 



The variables {u^,v^) represent "global" concentrations away from the perturbation and 
the concentration at the local perturbation. Tracking the growth or decay of this per- 
turbation provides stability information for ([T]). There arc two primary benefits to this 
approach: 1) the large amplitude "probe" detects pattern formation in linearly stable pa- 
rameter regimes, 2) the reduction to ODE's greatly simplifies its implementation. 

Applications of this method to biologically motivated reaction diffusion systems arc found 
in [131 [131 [12] • Rather than focus on a specific phenomena or biological system, my goal 
here is to explain the method itself. I will describe the types of RDE's to which this method 
is applicable, its limitations, and the type of information that it can provide. Well known 
examples of pattern forming systems are used to demonstrate its application and make direct 
comparisons between its predictions and results of classical methods (e.g. LSA, full PDE 
bifurcation, or numerical simulation) . In the context of a more complex chcmotaxis related 
example, I also show this method: easily scales to larger systems with many variables, allows 
the user to gain a more complete overview of the parameter space structure than with other 
methods, and greatly aids investigation of both parametric and structural perturbations of 
a complex reaction network. 
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2. Local Perturbation System Formulation 

I now proceed to show that the evolution of a spatially localized perturbation of a ho- 
mogeneous steady state of Eqn. (P) evolves according to Eqs. ©. Consider Eqn. ([T|) on the 
interval [—1,1] with no flux boundary conditions, and u,v in R^^ and M.^ respectively. It 
is not necessary to assume all slow (respectively fast) variables have the same diffusivities, 
only that they can be divided into fast and slow diffusing classes. For notation simplicity 
however, assume 

(3) Du = e^I, Dv=I?I, 

where I is the properly sized identity matrix. The central assumption will be that the three 
timescales defined by reaction kinetics, slow, and fast diffusion respectively arc substantially 
different, i.e. ^ 1 ^ D. Non-dimcnsionalization by domain size and the reaction 
timescale (so that f,g ^ 0(1)) have been implicitly assumed. Further assume this system 
has a HSS {u'^ {p) , v'^ {p)) satisfying f{u'^{p),v^{p);p) = = g{u'^{p),v^{p)\p). 
Consider a highly localized perturbation of this steady state of the form 

(4) ^t(a;,0)=^t^ ^;(a;,0)=^;^ |a;| > ^e, 
u{x,Q)=u^, w(x, 0)=w^, |a;| < ^/e, 

where {u^ ,v'p) is 0(1) with respect to e and D. Denote to be the local region \x\ < -y/e 
and the global region \x\ > -y/e. 

2.1. Time and space scale separation. To track the evolution of {u,v) on these regions, 
different time and space scales must be considered. A multiple timescale argument is applied 
to parse the role of reaction and diffusion effects on different timescales, and a boundary 
layer technique is used to separate relevant space scales. 

The reaction, slow, and fast diffusion timescales inherent in this class of RDE's can be 
described hy t = 0(1), = e^t, and ty = Dt, with the intermediate reaction timescale of 
most interest here. Suppose u = U{x,t,ty,ty),v = V{x,t,tu,ty). With the perturbation 
(HI taken as an initial condition, boundary layers on an 0(e) length scale are expected. 
Employing a stretched coordinate ^ = {x — xiayer)/^ near boundary layers, we come to two 
systems: 

(5) ^ = / + e'C/,„ ^=g + DV,, 

describing outer regions away from boundary layers and 

dU dV D 

(6) + 

describing dynamics in boundary layers. I now describe the evolution of (U, V) on the outer 
regions i?^'' on the short and intermediate timescales. Boundary layer effects and long 
timescale evolution are beyond the scope of this exposition. 

2.2. Evolution on the fast diffusion timescale. Consider first the fast diffusion timescale 
and assume U,V are described by first order pcrturbative expansions U = U'^ + eU^, 
V = V'^ + eV^. Substituting t^, = Dt into Eqn. ([S]) and collecting leading order terms 



(7) 
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it is clear that in outer regions, = U^{x,t,tu) does not evolve due to either reaction or 
diffusion and V'^ simply spreads due to diffusion 



(8) 



V°{x, t, tu, U) = v°{t, tu) + v"{t, tu) cxp("(ri7r)^t„) cos(ri.7rx). 



So in each outer region i?''^, evolves to a constant value exponentially quickly with t^. 



■ Outer Region 

■ Transition Layer 



0(Ve) 



-1 



0(e) 



Figure 1 . To probe non-linear stability, the Local Perturbation Analysis 
probes the response of a homogeneous steady state of Eqn. ([T]) to a localized 
perturbation Eqn. (jH). To leading order, values in the perturbed region and 
the broader domain, denoted it' , respectively, will evolve independently. 
The fast variable v, not depicted here, will be spatially uniform on the 
entire domain taking value . Boundary layer effects are present on 0{e) 
regions but to leading order do not influence evolution of the perturbation. 
A collection of ordinary differential equations for (u^,w^,u') describe the 
growth or decay of this perturbation and provide stability information for 

©■ 



Thus V'^) will evolve to a piecewise constant profile with possibly different values in 
and i?', see Figure [TJ Denote the values on the global domain by (u^, v^) and those 
on the local region i?' as (u',?;'). Now consider the boundary layers between these regions 
on this fast timescale. Given the spatial symmetry, consider only the left boundary layer 
and substitute ty = Dt into Eqn. 0{e~^) terms indicate that = 0. Thus to leading 
order 

^ao{ty)C + ai{t,). 

Matching conditions dictate that 

lim ^ v\ lim V"{0 = . 

So ao = 0, resulting in a shadow system |33l I26j where V*^ is constant over the entire 
domain. Its value will be denoted v^{t,tu). On this timescale to leading order, = so 
that boundary layer effects do not influence U°. 
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2.3. Evolution on the intermediate reaction timescale. Evolution of this perturba- 
tion on the reaction timescale will determine the stability of the HSS. As ty progresses 
approaching this timescale, to leading order the solution in the outer regions is described by 
{u^ ,u'-), representing the piecewise constant values on i?''^ respectively. The evolution 
of these values on this timescale are described by Eqn. ^ with t = 0{1). To leading order 

(9) ^^f{U°,V';p), ^^giU^,V';p). 

Substituting {u^,v^) and {u^v^) respectively into the first of these, we obtain Eqs. (|2ap . 
(Pc)) . Integrating the second of © over the domain we see that 



Qya I /-i 

(10) J ^9{U°,v<^;p)dx = g{u<^,v<^;p) + V~^[g{u\va;p) ^ g{ua,va-p)]. 

So to leading order the values (it',M^,w^) evolve according to ^ on the intermediate 
timescale. Eqn. ([2]) will be referred to as the LPA system of ODE's (or LPA-ODE's) as- 
sociated with Eqn. ([T]). Evolution on the slow diffusion timescale requires consideration of 
boundary layer effects. This is dependent on the specific system being investigated and is 
not discussed here. 

3. The Local Perturbation Analysis: Application and examples 

The goal of the LPA is to determine in which parameter regimes localized perturbations 
of this form grow or decay. Growth corresponds to instability of the HSS and a patterning 
response, decay corresponds to stability. Since the evolution of this perturbation is described 
by Eqn. ([2]), the location and stability of its steady states / fixed points describe stability 
properties of the HSS of Eqn. ([TJ. This information can be found using standard bifurcation 
analysis techniques for systems of ODE's. 

In coming sections, I will demonstrate this method through example and the following 
capabilities will be emphasized. 

(1) The LPA detects linear instabilities of Eqn. ([T|). 

(2) The LPA detects inherently non-linear patterning where a HSS is linearly stable but 
a sufficiently large perturbation yields a patterning response. 

(3) While the LPA approximation is not valid on the timescale of pattern evolution, its 
results can be used to make reasonable conjectures about the type of pattern (i.e. 
highly localized spike or a sharp interface separating distinct planer regions) that 
might evolve. 

(4) In non-linear patterning regimes, the LPA qualitatively maps the dependence of 
patterning response thresholds on system parameters p (excluding diffusion param- 
eters). 

The LPA is applied to two classical systems, Schnakenberg [35] and Substrate Inhibi- 
tion [22, both well studied in [31] |43l [T^ for example. Predictions of this analysis are 
then directly compared to results of linear stability, numerical, full PDE bifurcation, and 
asymptotic analyses. 

3.1. The local perturbation analysis of a Schnakenberg model. The Schnakenberg 
system is a Turing model where u is an activator and v a substrate. 

(11a) Ut {x,t) = a — u + u^v + e^Au = f{u, v) + e^Au, 

(lib) vtix, t)^h-u^v + DAv = g{u, v) + DAv. 
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u decays linearly, both are produced uniformly in the domain, and the nonlinearity represents 
an autocatalytic reaction where u consumes v. A LSA of Eqn. can be found in |31j . 
In [331 HI]: asymptotic and spectral analysis showed that for a = 0, highly localized spike 
solutions exist and arc stable. 

It is possible to analytically perform the LPA for Eqn. pT|) . The resulting system of 
LPA-ODE's becomes 

(12a) ul^a-u3 + {uafv\ 

(12b) v1 ^b^{u3fv\ 

(12c) u\^a-u^ + {u^fv3. 

Here, p = (a, h) is the vector of system parameters and a will be the bifurcation parameter 
of interest. Eqs. (|12ap . (jl2bp decouple from Eqn. (|12cp and simply represent the spatially 
homogeneous system (i.e. with e = = D). The unique HSS [u'^ [p) , [p)) of Eqn. (|lip 

(13) u^=a + b, 1'^ = 7 — TTT^' 

(a + b)'^ 

is thus a solution of Eqs. (|12ap . (jl2bp . Similarly {u^,vS,u'') = (u*,w'*,w*) is a steady state 
of Eqn. (HI]). This steady state of the LPA-ODE's represents the HSS of Eqn. ^ with no 
perturbation, i.e. both u'-^ = u'^. 

While this is the only HSS of the RDE system, the LPA system actually has two steady 
states with the second satisfying 

2 

(14) u3 = u\ = v\ = a + — u'^ 



With b fixed and a considered as a bifurcation parameter, the steady state branches and 
m'^ intersect in a transcritical bifurcation at a = &. Furthermore, it can be readily shown 
by computing the Jacobian of Eqn. that the stability of these branches is determined 
solely by the sign of /„ and that 

(15) ^iu',vn>0, ^{u'\vn<0, a<b, 

(16) ^{u',vn<0, ^{u'\v')>0, a>b. 

The location and stability of these steady states is depicted in Figure [2^. The HSS branch 
{u^,vS,u'') = {u'',v^,u'') is linearly unstable for a < 1 (Region I). For a > 1 (Region II), 
the HSS is linearly stable, however a perturbation of above the u'^ branch will grow to 
infinity. From here on, the HSS branch {u^,v^,u'') = {u^,v^,u^) will be referred to as a 
"global" steady state branch of the LPA system. {u^,v^,u'') = (it'', w*, it'^) will be referred 
to as a "local" branch, since it is only a steady state for the local variable uK 

3.1.1. Local perturbation analysis predictions. These results lead to the following predic- 
tions. 

Prediction 1: A Turing bifurcation occurs near a = 1. For a < 1, the homogeneous 
steady state of Eqn. pT|) is linearly unstable. For a > 1, it is stable but sufficiently large 
perturbations yield a patterning response. 

Based on the asymptotics above, it is expected that the initial behaviour of a local per- 
turbation of the RDE's mimics the behaviour of the perturbation u' determined by this 
bifurcation analysis. In region I (a < 1 in Figure [2^), arbitrarily small perturbations of 
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= grow, predicting the HSS of Eqn. ([TT|) is linearly unstable. In region II, sufficiently 
large perturbations of are required to elicit a response for the LPA-ODE's, suggesting 
the HSS is linearly stable but large perturbations yield a response. 

Prediction 2: In region II, as a increases, increasingly large perturbations are required to 
initiate patterning. 

In region II (a > 1), the gap between the stable global and unstable local branches represents 
a response threshold for the LPA-ODE's: a perturbation of below the threshold decays 
back to the global HSS branch, a perturbation above it grows. The dependence of this 
response threshold on the system parameter a can be found by visual inspection of the LPA 
diagram. For a > 1, that threshold increases with a. This threshold is precise only in the 
e — > 0,D — >■ oo limit, however the qualitative dependence on a is expected to hold for the 
RDE system pTjl with sufficiently extreme diffusivities. 

Prediction 3: The Turing bifurcation near a = 1 is sub- critical for sufficiently extreme 
diffusivities with large amplitude patterned states present on both sides of the bifurcation. 

In dynamical systems theory, the terms sub-critical and super- critical are often used to 
describe the character of bifurcations such as Hopf or pitchfork. Super-critical denotes 
a bifurcation that gives rise to a small amplitude response upon crossing it. Sub-critical 
denotes one where the HSS loses stability immediately giving way to a large amplitude 
response. In the latter case, responses can occur even outside of the unstable regime given a 
sufficient perturbation. The Turing bifurcation near a = 1 is predicted to be sub-critical with 
the unstable local branch u'^ characterizing a threshold that shrinks to at the bifurcation, 
giving rise to instability of the HSS. 

Prediction 4: Patterned solutions take the form of a spatially localized spike. 

In region II, the LPA-ODE's exhibit blow up; a perturbation above the critical threshold 
grows to infinity. When this perturbation becomes large, diffusion is expected to become 
important for the RDE's. This will tend to oppose growth and smooth the solution. It is 
reasonable to conjecture that at a particular height, reaction driven growth and diffusion 
driven suppression will balance leading to a large amplitude, spatially localized spike. For 
future reference, it is expected that when e decreases, decreasing the strength of diffusion, 
this spike would be come taller and more localized. 

3.1.2. Verification of predictions. Figures [2] (b,c,d) show results of linear stability, full PDF 
bifurcation, and numerical analyses for Eqn. (jlip . These results support the predictions 
above with a few caveats discussed at the end of this section. 

Verification of prediction 1 

Detection of a Turing instability for a < 1 agrees well with LSA results presented in Figure 
[2)d. There, eigenvalues of the linearized Jacobian (i.e. Turing growth rates) are plotted as 
a function of the bifurcation parameter a for four successively smaller values of e. Dots on 
Figure [2b indicate LSA bifurcations; these bifurcation values are recorded for two different 
values of D in Table 13.1.21 The location of these bifurcation values appears to converge to 
the predicted value of a = 1. This supports point 1 in Section |3] that the LPA detects linear 
instabilities. 

The LPA predicts that in the linearly stable region II, sufficiently large perturbations elicit 
a patterning response. To test this, both full PDF bifurcation analysis and asymptotics are 
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Schnakenberg LPA Schnakenberg LSA 




(c) (d) 

Figure 2. Comparison of linear stability, local perturbation, and full PDE 
bifurcation analysis results for the Schnakenberg system (|lip . All diagrams 
are computed with D = 10, b = 1. Panel a) Local perturbation anal- 
ysis results: Global and local u'^ solution branches of Eqn. along 
with their stability are plotted as a function of a. Two pattern forming 
regimes are predicted; I) linearly (Turing) unstable, and II) linearly sta- 
ble where a sufficiently large perturbation induces patterning. Panel b) 
Linear stability analysis results Maximum eigenvalue of Ji (fT8|) as a 
function of "a" for various values of e. As e 0, the edge of the Tur- 
ing region, marked with dots, approaches a limiting point near a ~ 1, in 
agreement with LPA predictions in panel a. Panel c) PDE bifurcation 
results: Bifurcation analysis of the full system of PDEs ([TT|). The vertical 
axis describes the height (maximum - minimum) of a patterned solution. 
As predicted by the LPA, stable patterned solutions exist both in the lin- 
early unstable and stable regimes for sufficiently small e. The location of 
the Turing bifurcation near a = 1 agrees with panels a,b. Marked points 
represent points where the computed solution is plotted in panel d. Panel 
d) Example solutions of Eqn. with e ~ .025. 

employed. Figure [2t shows results of numerical continuation (using Auto [6]) of patterned 
solutions of the full PDE system (jlip with the vertical axis depicting the height (maximum 
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Table 1 . Value of "a" at the edge of the Turing region for the Schnaken- 
berg system (jlip with 6=1 for various values of e. Column two: values 
drawn from the marked points in Figure [21d. Column three: similar values 
for D = 1000. The final row is the value of the bifurcation predicted by 
the LPA. This bifurcation approaches that predicted by the LPA as e — ^ 0, 
13 ^ oo in agreement with Corollary 14.21 



e 


Turing (D = 10) 


Turing (D = 1000) 


0.1 


0.76 


0.82 


0.05 


0.88 


0.95 


0.025 


0.91 


0.98 


0.01 


0.93 


0.99 


LPA 


1 


1 



- minimum) of the patterned solution. The horizontal axis depicts the unpatterned HSS 
(maximum - minimum=0). For sufficiently small values of e, the stable patterned solution 
extends into the a > 1 region where the HSS is linearly stable, verifying the presence of stable 
patterned solutions in that regime. Further, asymptotic results in Appendix |A] show that 
in the e — >■ 0, 13 — >■ oo limit, patterned solutions exist for all values of a. This lends support 
for point 2 in Section |3] that the LPA detects inherently non-linear patterning regimes. 

Verification of prediction 2 

In Figure|3]local perturbations of different height were applied to the HSS for multiple values 
of a and the presence / absence of a pattern was recorded. As a increases, the perturbation 
size required to induce patterning increases as predicted by the LPA. This supports point 
4 in Section [3] that results of the LPA can be used to determine the qualitative dependence 
of response thresholds on parameters in non-linear patterning regimes. 

Verification of prediction 3 

Figure[2j: shows that as e — )■ 0, the nature of the Turing bifurcation changes from being super- 
critical to sub-critical. For large e, small amplitude patterns emanate from the bifurcation. 
For smaller e, the stable HSS gives way to large amplitude pattens immediately upon crossing 
the bifurcation. Also for small e, an unstable patterned state is present outside the linearly 
unstable regime. As the bifurcation is approached, this unstable state collapses onto the 
HSS changing its stability. This is similar to the standard example of a sub-critical Hopf 
bifurcation where an unstable limit cycle colliding with a stable node results in a bifurcation. 

It seems this association of a non-linear patterning regime with a sub-critical Turing bi- 
furcation is somewhat common. Both the substrate inhibition example and the chemotaxis 
example in Section [5] presented later show sub-critical bifurcations. Similarly, Rodrigues et 
al. [35] observed the presence of a stable heterogenous state adjacent to Turing parameter 
regimes for a discrete predator prey model. They also showed the response threshold neces- 
sary to induce patterning in these regimes decreases as the Turing bifurcation is approached, 
consistent with these results. 

Verification of prediction 4 

Figure [2ji and results in Appendix |X] show that in both the linearly stable and unstable 
regimes a spike like solution forms. Furthermore, these results show that as e is decreased 
(reducing the opposing effect of diffusion) , the spike height increases as expected. Thus the 
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Figure 3. Verification of tlie qualitative relationsliip between "a" and tlie 
patterning tliresliold in tlie Schnakenberg system using numerical sim- 
ulation, e = .01 and all other parameters are as in Figure [21 Simula- 
tions were given a period of time to settle into a stable steady state (when 
present). Perturbations of varying size were then applied to the middle 10% 
of the domain. The vertical axis represents the size of the applied perturba- 
tion, 'x' indicates the perturbation grows resulting in a spike, 'o' indicates 
the perturbation decays back to the homogeneous state. The patterning 
threshold increases with a as indicated by the LPA. To the left awl, the 
homogeneous state is unstable and to the right of a « 1.7, it is stable to all 
perturbations, in agreement with Figure [2j;. 



LPA inferences about the long term evolution of patterns are confirmed in this example, 
supporting point 3 in Section [S 



3.1.3. Caveats of the local perturbation results. There are discrepancies between the LPA 
predictions and results of LSA and full PDE bifurcation analyses. First, the location of the 
predicted bifurcation at a = 1 is not precise. Both LSA and PDE bifurcation results show 
the value of the actual Turing bifurcation depends on e and D. Though this does appear to 
converge to the predicted a = 1 in the proper limit. This type of approximation error will 
be present in any LPA application and will be discussed in more detail in Sectional 

Second, the LPA predicts patterned solutions will form for all a > 1. Full PDE bifurcation 
results (Figure [2}:) in contrast show the patterned state is annihilated in a saddle node (or 
fold) bifurcation at a finite value of a = a*{b,e,D). The location of this bifurcation does 
increase as e — >■ and asymptotic results in Appendix |A] show the presence of a patterned 
solution for all values of a. So in the e — > 0, D — > oo limit, a* — > oo, in agreement with LPA 
results. 



3.2. The local perturbation analysis of a substrate inhibition model. I now ap- 
ply the LPA to a substrate inhibition model to demonstrate a different set of results and 
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interpretations obtained with the same method. This model [21] 

(17a) ut(x,t) = a — u ^ + e^Au = f(u, v) + e^Au, 

1 + u + Kw^ 

(17b) vt{x, t) = a{b -v) , + DAv = g(u, v) + DAv, 

1 + u + A v/ 

describes two co-substrates that are constantly generated, decay linearly, and are used up 
in an enzymatic reaction. The non-linear term is indicative of multiple substrate molecules 
u binding to a single enzyme rendering it inert for further interaction with the remaining 
co-substrate v, thus the term substrate inhibition. See [31] for LSA results for this system. 

In the previous example, it was possible to analytically compute the various solutions of 
the LPA system of ODE's. This will not generally be the case, but it is possible to find 
and track the various LPA solution branches efficiently using standard ODE bifurcation 
techniques. Figure |4] mirrors Figure [2] with results in panel b,c,d verifying predictions in- 
ferred from LPA results in panel a. In this example, solution branches of the LPA-ODE's 
along with there stability are computed with the numerical continuation software pack- 
age Matcont LPA predictions regarding the presence of linearly unstable (region III) 
and non-linear (regions II, IV) pattering regimes are again verified by LSA and full PDF 
bifurcation analyses. 

There are a few important contrasts between these LPA results and those for the Schnakcn- 
berg model. First, in region IV (Figure S^), a negative valued perturbation of the HSS is 
predicted to induce patterning. This along with the general dependence of response thresh- 
olds in regions II, IV were verified numerically (results not presented). Second, the LPA 
results suggest the resulting solution will take the form of a stable interface separating high 
and low regions for u. 

Consider region II where this system has two local branches, one unstable and the other 
stable. In this case, a perturbation of above the unstable local branch will be attracted 
to the stable local branch. On the slow diffusion timescale, boundary layers between the 
high and low u states on respectively will move. The lack of a higher HSS suggests 
the high state cannot encompass the entire domain and one of two things will happen: 1) 
the boundary layers will stop / stall somewhere in the interior of the domain leaving stable 
interfaces separating regions of high and low u levels, or 2) the solution will eventually 
collapse back to the unique HSS. Figure [4(d)] showing stable long term solutions for multiple 
values of a verifies the former is the case here. These results again support the suppositions 
in Section [U 



4. Detection of linear instabilities by the LPA 

The previous section demonstrated points 1-4 in Section [3] Points 2-4 relate to non-linear 
patterning regimes and will be difficult to show in generality. The supposition in point 1 does 
however hold in generality. In this section I show the LPA detects linear instabilities with 
increasing accuracy as diffusivities become extreme. Asymptotic analysis shows eigenvalues 
of the linearized RDE's converge to eigenvalues of the linearized LPA-ODF's. Thus positive 
eigenvalues of the LPA system (indicating linear instability) correspond to positive Turing 
growth rates and vica versa. 

Since eigenvalues of the linearized RDE's ([Ij determine linear stability properties, they 
will be characterized in the e — >■ 0, -D — >■ oo limit. Linearizing about the HSS (u^,w^) with 
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Figure 4. Comparison of linear stability, local perturbation, and 
full PDE bifurcation analysis results for the substrate inhibition 
model Eqn. p7|) . All diagrams are computed with D = 10, p = 13, 
K = 0.125; a — 1.5, and 6 = 80 and all conventions are as in Figure 
[21 Panel a) Local perturbation analysis results: Three regimes of 
behaviour are predicted: No patterning (1), linearly (Turing) unstable (111), 
and linearly stable where a sufficiently large perturbation yields a response 
(II, IV). Panel b) Linear stability analysis results: For each e two 
Turing bifurcations arc seen. Dots mark the location of the right Turing 
bifurcation for different values of e. As e — >■ 0, the edges of the Turing region 
approach that predicted by the local perturbation analysis results in panel 
a. Panel c) PDE Bifurcation Results: In agreement with the LPA 
results, patterned solutions are found both inside and outside the linearly 
unstable regime. Panel d) Example solutions drawn from starred points 
for e ~ .05 in panel c. As predicted by the LPA results, these solutions 
show a stable interface separating high / low regions of u. 



respect to periodic perturbations of the form cxp(iA:a;) leads to the Jacobain 
/1o^ r/„(u",u'';p) - Pe^/ /„ (u" , u'' ; p) 
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Recall that / : R*^ x ^ R^^ , g : R^^ x R^ ^ R^ and denote eigenvalues of this matrix 
as {A^(e, i^,p)}i=i:(A/+Af)- Assmiie these are in decreasing order according to their real part 
so that \\ has the largest real part and determines stability. 

Stability of the HSS branch of the LPA ODE's (u^ , , v}) = (it'*, w", u'*) is determined by 
eigenvalues of the LPA Jacobain 



(19) 



Jlp 



fu{u^,v^;p) 



This is a block lower triangular matrix with the upper left block being precisely Jg (|18p . 
Thus the eigenvalues {A°}i of Jq are also eigenvalues of Jlp- These describe stability of the 
HSS of Eqn. ([IJ with respect to spatially homogeneous perturbations and determine stability 
of the well mixed system (i.e. with I?„ = = Dy). So the LPA detects any instabilities of 
the well mixed system. These will be referred to as the well mixed eigenvalues of Jlp- 

The remaining eigenvalues of Jlp are simply eigenvalues of /„(w^,w^). Denote these as 
{^j'^{p)}j=i:M and assume they are ordered according to decreasing real part so that Af ^ 
has the largest real part and determines stability. Then the following asymptotic result 
relating {A^^} and {A'^} holds. 



Theorem 4.1. Assume ^ 1 <^ D, V/, Wg are 0{1) with respect to e and D, and fix a 

wave number fc > 0. Further assume that /^(w'', v'*), 5t,(u'*, v'*) are diagonalizable. Then: 

(1) For each i = 1 : Af = Af ^ - fc^e^ + c{D) where c{D) as D ^ oo. 

(2) For each i = AI + I : AI + N , fHe(Af) — 0^{D) where 0^{) signifies a negative 
valued quantity of that order. 



For proof of this result, sec Appendix |B] A direct consequence of this is that as e 0, ^■ 
oo, Xi{e,D,p) Xi^{p) independent of fc. So in this limit, linear instability of the HSS 
branch of the LPA-ODE's corresponds directly to diffusion driven (i.e. Turing) instability 
for the RDE's. Thus point 1 in Section [3] holds for general systems of the form ([T|) when 
fu{u'',v'^),gy{u'',v'') are diagonalizable. 



4.1. LPA bifurcations locate the edge of "Umiting" linearly unstable parameter 
regimes. Recall from the LSA results in the previous examples that the location of a Turing 
bifurcation of the RDE's appears to converge to a limiting point as e ^ 0, £> — ^ od. This 
is in line with Murray's |31j observation that a linearly unstable regime of parameter space 
converges to a "limiting" unstable regime in this limit. As indicated by the comparison of 
the locations of LSA bifurcations and those predicted by the LPA, the LPA precisely locates 
the edge of these "limiting" unstable regimes. This follows from the following corollary of 
Theorem O 



Corollary 4.2. Consider a particular set of parameters p. If the global branch of the 
LPA-ODE's is linearly unstable ( i.e. fRe(Af^(p)) > Q), then the HSS of the RDE's is 
linearly unstable for sufficiently extreme values of e and D (i.e. lHe(A^(e, I?,p)) > for 
some k > 0). Furthermore, if the global branch is linearly stable in the LPA sense, the HSS 
is linearly stable for sufficiently extreme values of e, D as well. 
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5. Applications of the LPA to more complex systems 

One of the primary benefits of the Local Perturbation Analysis is the relative ease with 
which it can be applied to more complex systems, common in biological applications. Ex- 
isting methods become either difficult to implement or intractable in such cases. However, 
with the help of ODE analysis software packages such as Auto [5] and Matcont [S], this 
method scales well to larger systems. The following example demonstrates an application 
of the LPA to a system involving 9 RDE's. 

5.1. Chemotactic polarization example. Much effort has been devoted to understand- 
ing the process by which cells, ranging from white blood cells to cancer cells, move up 
chemical gradients. Reorganization of regulatory molecules, primarily GTPases and phos- 
phoinositides, is known to be a precursor to such motion. In response to an applied chemical 
gradient, these molecules self organize to form a polar state where some localize in the cell 
"front" (Cdc42, Rac, PI3K, and PIP3) and others in the "rear" (Rho, PTEN). Front re- 
lated molecules generate protrusion, rear related molecules generate contraction, and their 
combined activity leads to directed motion. 

Each of the three GTPases (Cdc42 C, Rac R, and Rho p) effectively has two forms, 
membrane bound and cytosolic with only the membrane bound form in an active state. 
Over the timescale of polarization events, the amount of each GTPase is conserved with 
diffusion and cross talk mediated cycling between the two states leading to segregation of 
active forms. These cross talk interactions and the influence of phosphoinositide feedback 
are the focus of this discussion. 

While these regulators are conserved across a wide range of eukaryotic cells, the cross 
talk interactions between them is not. This variation has led to extensive experimental work 
aimed at dissecting these interactions in different cell types and numerous models (reviewed 
in |19j ) aimed at understanding their results. Here I describe and analyze a variant of a 
model [UJ [37] motivated by work on HcLa cell polarization. 

I investigate a structural perturbation of that model, introducing mutual antagonism 
between Rac and Rho, known to be present in numerous cell types [37l[2l|42]. A schematic 
diagram of this model is in Figure 15.2b .. The dashed interaction, Rho mediated inhibition 
of Rac, is the structural addition differentiating this model from that in [Ml [27]. Model 
equations encoding these interactions are found in Eqs. (|34p with a description of parameters 
and their values in Table [5] See Appendix [C] for a brief description of this model and [Hl[?7] 
for more extensive discussion of the original network and its parameters. 

What effect does the addition of Rho mediated inhibition of Rac have on the behavior 
of this network? To investigate this, a non-dimensional parameter (/2) modulating the 
strength of this inhibitory interaction is introduced. When /2 = 0, no inhibition is present 
and the original network is recovered. When /2 increases, the strength of the inhibition 
increases. LPA and numerical simulation results in Figure show the effect of increasing 
the strength of this feedback. 

5.2. LPA and numerical simulation results. Figure 15.2b shows the results of a LPA 
of this model with moderate feedback, /2 = 2. The LPA was performed assuming mem- 
brane bound GTPases are slow diffusing = O.lfim^ / s), and cytosolic GTPases {Dc — 
BOfim^/s) are fast. For reference, cell sizes considered are on the order of 10 — 20/im. Phos- 
phoinositide diffusion lies between the fast and slow regimes. However, as in jl4j . LPA 
results are similar with it chosen as either fast or slow. In Figure 15.21 they are taken to be 
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slow variables. The LPA reduction in this case leads to a system of 15 ODE's for 6 local 
variables {C\R\p\Pl P^, Pj) and 9 global variables {Cs.Rs, p9,Ca,R9, pl.Pf, Pi, Pi)- 

The bifurcation parameter Ijn represents a basal Rac activation rate. Its variation could 
result from either population heterogeneity or external stimulation of Rac as in [57]. In 
Figure 15.2b . three regimes of behavior are found at different activation levels. For both 
low and high levels of basal activation, no response due to either instability or an applied 
stimulus can occur. For increasing levels of activation, a regime where sufficiently large 
perturbations yield a response is found. Again, this regime terminates in a sub-critical 
Turing bifurcation as the response threshold shrinks to zero. This suggests increasing the 
spatially uniform activation rate increases the sensitivity of a cell to heterogeneous stimuli, 
experimentally supported in [?7] . 

At yet higher values of Im, a second narrow linearly stable patterning regime is found. 
Numerical simulations verify the presence of all but this narrow regime, which likely requires 
more extreme diffusivities to be observed numerically. LPA results again suggest solutions 
will evolve to a polarized profile with a transition layer separating regions of homogeneous 
activity levels. This was also verified numerically with an indicative steady state solution 
shown in the inset (Iri — 1.1). 

Now consider the effect of increasing / decreasing the strength of Rho H Rac feedback. 
Labeled branch points (BP) indicate the approximate boarder of the unstable regime. Fold 
bifurcations (LP) mark the boarder of the non-linear patterning regimes. Standard two 
parameter continuation techniques are applied to follow these bifurcations as /2 is varied, 
Figure [Ob . At low values, both regimes persist. As /2 increases, the branch points collapse 
at /2 ^ 5 and the linearly unstable regime between them is lost. For higher values of /2, 
the two fold bifurcations of the local branch persist suggesting the continued presence of a 
linearly stable patterning regime between them. 

Marked points on Figure[52t indicate parameter values where numerical simulation of the 
full RDE system was performed. Circles indicate a parameter set where small noise (machine 
noise) induces a response. Points marked x indicate sufficiently large perturbations are 
required for a response. Beyond these points, no patterning was detected numerically at 
the base diffusion values. When Dm = .01/xm^/sec, parameter regimes expand with squares 
marking additional parameter sets where sufficiently large perturbations yield a response. 

The linearly unstable regime for the RDE's is confined to that predicted by the LPA 
and as expected, stimulus induced patterning is present to the left but not to the right of 
that regime. While the location of patterning regimes in parameter space agree well with 
predictions, the expanse of these regimes is substantially smaller than predicted, particularly 
for the non-linear patterning regime. However, as diffusivities are driven to yet further 
extremes, these regimes do expand further (results not presented). 

I consider one final perturbation of this network, PI3K knockout which is accomplished 
here by setting kpj3K=o- This has the effect of removing feedback between GTPase's and 
phosphoinositide's. The same analysis above was performed with results in Figure [5T2l l. Sim- 
ilar parameter space structure exists with linearly unstable and stable patterning regimes. 
In this case however they are substantially compressed in parameter space. So while PI3K 
/ PIPS mediated feedback is not necessary for polarization, it does make it more robust 
in a parametric sense. This is consistent with observations [H [27] showing PI3K / PIPS 
localization is not necessary for efficient chemotaxis, but its knockout substantially reduces 
the fraction of cells that do clicmotax. 
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(c) (d) 

Figure 5. Panel a Schematic diagram of interactions between three GT- 
Pases (Cdc42, Rac, and Rho) along with phosphoinositides (PIPl, PIP2, 
PIPS) represented by Eqs. (p4| . An arrow represents activation of one pro- 
tein by another, a bar represents inactivation. Panel b: Local pcrtmbation 
analysis of this system with ji = 1 and l^x the bifurcation parameter. The 
vertical axis is the value of the local form i?' of Rac. The monotonic branch 
is the global branch, the loop is the local branch Four regions arc seen, or- 
dered from left to right: stable, linearly stable where a sufficiently large 
perturbation yields a response, linearly unstable, and stable. There is a 
small nonlinear patterning regime regime between l^x ^ 1.6,1.7, but this 
regime appears to only be present for extreme diffusions beyond those con- 
sidered here. Inset: Numerical simulation results at ji = 2, Ij^i = 1.1. 
Panel c: Two parameter continuation of the branch points (BP), indicat- 
ing the edge of a linearly unstable regime, and fold bifurcations (LP) of the 
local branch. Markers represent simulation results of the full RDF system. 
Circles indicate a Turing instability, 'x' indicates perturbation induced pat- 
terning, and a square indicates a parameter set which is stable for D„i ~ .1 
but where perturbation driven patterning results with Dm = -01. Panel d: 
The same as Panel c with PI3K knockdown removing the feedback of PIPS 
— > Rac. Similar parameter space structure is present though substantially 
compressed. Values for parameters not explicitly mentioned are in Table [5) 
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6. Limitations of the LPA approximation 

Here I stress the limitations of the Local Perturbation Analysis. The purpose of the LPA 
is not to approximate a solution of Eqn. ([T|), only the initial response of a HSS to a localized 
perturbation i.e. growth or decay. Slow diffusion timcscale and boundary layer effects 
are not considered and the LPA-ODE's ^ only describe the evolution of the perturbation 
on short to intermediate timescales. These effects can become important and the leading 
order approximation above can fail for one of two basic reasons. 

• The perturbation becomes large. This would cause a number of effects. First, if g is 
unbounded, the correction term in Eqn. (fTO| could become 0{1) and affect leading 
order dynamics. Second, if /, g arc unbounded, neglected Taylor expansion terms 
of the form e/„, e/„, eg„, or eg^ can become 0(1), influencing the leading order 
dynamics. Third, boundary layer effects could influence the dynamics of [U, V) on 
RK 

• The perturbation spreads in space. This would again cause the area of the perturbed 
region to become 0(1), causing the correction term in Eqn. pUj) to affect leading 
order dynamics. These effects would however occur on the slow diffusion timcscale 
which is not considered here. 

In either case, these effects only become important after the perturbation has grown in 
amplitude, constituting a response. 

As a result of neglecting these higher order effects, LPA predictions are asymptotic in 
nature and require sufhciently distinct diffusivities to be valid. Therefore the predicted 
location of bifurcations between parameter regimes only approximate the location of actual 
bifurcations, with this approximation improving as diffusivities become more extreme. In 
some cases, bifurcations of the RDE system with finite diffusivities (such as the saddle node 
bifurcation in the Schnakcnbcrg case) are not captured at all by the LPA. 

Further, all diffusion related information is lost. Therefore length scale information can- 
not be obtained, non-linear phenomena such as peak splitting [531 11] will not be found, 
and pattern selection or abbarent effects from domain growth for example [H [T] cannot be 
discussed. For these reasons, the LPA should be viewed primarily as an efficient, scalable 
non-linear stability technique that can be used to inform further analysis. 

7. Discussion 

A new non-linear bifurcation technique for systems of reaction diffusion equations with 
large diffusion disparities ^ was developed and demonstrated. This "Local Perturbation 
Analysis" (LPA) determines the response of a HSS of a system of reaction diffusion equations 
to a spatially localized, large amplitude perturbation. Under proper asymptotic assump- 
tions about the diffusivities e, D and the form of the perturbation, its evolution can be 
approximated to leading order by a collection of ODE's describing the perturbation (local 
variables) and the broader domain (global variables). 

A bifurcation analysis of this collection of LPA-ODE's reveals two types of solution 
branches : 1) a "global" branch of solutions representing HSS solutions of the RDE's, 
and 2) "local" solution branches unique to the LPA-ODE's. The location and stability of 
the global branches provides linear stability information for the RDE's ([T]). The location 
and stability of the local branches provides non-linear stability information. Application of 
this method and interpretation of its results were demonstrated using two classical example 
systems, Schnakenberg and substrate inhibition. To demonstrate scalability of the LPA to 
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larger more complex systems, it was applied to a biologicaly motivated system involving 
nine interacting chemotaxis regulators. 

The Local Perturbation Analysis has a number of benefits, i) It is simpler to implement 
than most PDE analysis techniques. These are notoriously challenging and typically tai- 
lored to the particular system being investigated, ii) With the help of existing software 
it can be applied to systems with many variables, common in biological applications, iii) 
In this setting, effects of both parameter and structural variations of a reaction network, 
resulting from cells using the same conserved biochemical machinery in different ways, can 
be efhcicntly investigated, iv) Since the method requires only analysis of ODE's, it should 
be accessible to a broad base of uses not versed in non-linear PDE techniques. 

Beyond these practical considerations, the LPA concisely summarizes a rich set of linear 
and non-linear information on a single one or two parameter bifurcation diagram, v) It accu- 
rately detects the location of linear instabilities (when diffusivities are sufficiently different), 
vi) It detects non-linear patterning regimes which appear to commonly be associated with 
sub-critical Turing type bifurcations, vii) In these regimes, it qualitatively characterizes 
the dependence of response thresholds on reaction parameters, viii) The global bifurcation 
structure can be interpreted to provide reasonable hypotheses about the type of pattern 
that will evolve on longer timescales. For these reasons, the LPA has the potential to be of 
use in an array of scientific fields where RDE's arise. 
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Appendix A. Schnakenberg asymptotics 

This analysis closely follows [33] . Consider the Schnakenberg system pT|) on the interval 
[— 1, 1]. In [43] it was shown that this system exhibits stable spike solutions when a = 0. 
That analysis can be extended to show such spikes in fact exist for all values of a under 
certain asymptotic conditions. 

To begin, define 

D _ u 

(20) D = — , V = ev, u = —, 

e e 

and subsequently drop the "to yield 

(21) Ut{x^ t) = ae — u + u^v + eu^^ 

2 

(22) evt{x,t)=b- — +Dv,,., 

e 

Assuming D ^ 1/e, v — vq + evi{x) + .... and integrating (j22p . it can be determined that 

(23) 26 = — / u^{x,t)dx. 



e 



1-1 

A spike solution of the form u{x) = ua + ui{x/e) is now sought where wq and ui are the 
outer and inner solutions. It is assumed uo is spatially constant. Collecting terms involving 
the same powers of e shows the outer solution is uq ~ ae and the inner solution satisfies 

(24) u'^{z)-ui{z)+ul{z)vo:^0 
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on xjt = z S [—1,1] with no flux boundary conditions. The solution to this problem is 
known (see |43| ) yielding 

3 X 

(25) u(a;) = uo + Ml = ae H sech^( — ). 

2wo 2e 

Integrating the square of this expression and substituting into (j23p yields vo = ^/3- Un- 
ravelling the change of coordinates yields the approximate spike solution for the original 
problem ((TTl) on [-1, 1] 

(26) u{x) w a + ^ scch^(^), v{x) w 

2e 2e o 

So the Schnakenberg system ([TT|) in fact produces spike type solutions for all values of a 
in the limit e — >■ 0. This is in agreement with the results of the LPA in Figure [5^ and the 
progression of the fold bifurcation (where the spike is lost) to cxd as a — oo in the bifurcation 
analysis in Figure [5};. Further, the maximum value of u in (j26p with a = compares to good 
precision with the maximum values shown at a = in Figure [2]:, supporting these results. 

Appendix B. Proof of Theorem 14.11 
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Figure 6. Schematic of the separation of the eigenvalues of Jfe in the 
complex plane. Grey circles indicate the different Gershgorin circles Ci. 
The larger darker circles indicate C; and Cg which contain all eigenvalues 
of Jfc. These circles separate those eigenvalues into two classes with 0(1) 
and 0~ (D) real part respectively. 

To prove Theorem 14. 11 first notice the eigenvalues of Jk ([T^ can be segregated into two 
regions of the complex plane using the Gershgorin circle theorem (see Figure [6]) . Fix a 
specific wave number k, let a^.j be the elements of Jk, and define 

(27) i?, =^|a,,,|, a = C{au,R,) 

where C{a,r) is the circle with centre a and radius r. The Gershgorin circle theorem states 
that each eigenvalue of Jk lies in at least one of the disks Ci . The structure of Jj, is such that 
the off diagonal entries are 0(1) with respect to -D. So i? = max{i?i} is 0(1). The diagonal 
entries fall into two categories, those that are 0(1) (corresponding to the small diffusion 
entries), and those that are —k^D + 0{l) (corresponding to large diffusion entries). Define 
Jls to be the union of the disks d that are characterized by 0(1) diagonal entries and ili as 
the union of disks characterized by 0{D) diagonal entries. Since these disks have a maximal 
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radius R independent of D, there exists disks C; = C{—k^D, kiR) and Cg = C{0, KsR) so 
that for constants k;_s independent of D, fi^ C Cg and fi; C Cj. For sufficiently large D, Cg 
and Ci do not overlap and hence separate {A*^} into two sets (see Figure IH]). 

So, for each i, either ^He(Af) = 0-{D), or ^HeCA*'') = 0(1). Since det(Jfc) = ©(Z?^), 
{A*^}i=A/+i:A/+Ar rnust have 0~ (D) real part. Also note that the imaginary parts of all 



eigenvalues are constrained to be less than max{Ks, ki}R so that Ofm(A^) 
well, so |Af I = 0(1) for i = l:M. 

Eigenvalues of Jk are roots of the characteristic polynomial 



(28) 



\Jk-Xl\ = 



A)/ Mu',v') 

9v{u',v') - {k^D + \)I 



0{1) for all i as 



= 0, 



where / is a properly sized identity matrix. Let P and Q be the unitary matrices that diag- 
onalize /„(u*,t;*) and gv{u'^,v^). Then in particular the diagonal entries of P^^/„(u'% tr')P 
are {Aj"^} and the entries of Q~^gv{u^ ,v^)Q are 0(1). Define 



(29) 



T = 



P 
Q 



Then the eigenvalue problem translates to 
(30) 



[A^^] - k\^I - XI 


P-'fvQ 




Al 


A2 


Q-^guP ( 


Q-^gv{u'',v'')Q- k^DI - \I 




A'S 


AA 







where [A^^] is the diagonal form of fu{u^,v'^). Notice that A1,AA arc diagonal. 

Now consider an eigenvalue A whose real part is 0(1). In this case, the diagonal entries 
of A4 arc 0~{D) and it is non-singular. It can thus be used to eliminate A2. After this is 
done, the eigenvalue problem becomes 



(31) 



[A 



k^e'^I + 0{D- 

Q-^guP 



XI 







Q-^gy{u',v')Q^k^DI 
Since the bottom right block is non-singluar, it must be true that 



XI 







(32) 



dct ([A 



LPi 



k\^I 



0{D- 



XI) 



0. 



where 0{D~^) is a properly sized matrix with entries of this size. With D = cx), the roots of 
this polynomial are simply {A^''^ — k^e^}. It is tempting to view Eqn. p2p as a perturbation 
of this case and apply some form of perturbation bound. However, /„ is not Hermitian, 
which is usually required for such bounds. Instead, the best we can say is that by continuity 
of the determinant, the roots of this polynomial satisfy 



(33) 

where c{D) 



A = Ar - k' 



as _D — )• oo. 



Appendix C. GTPase model equations 

Figure inS^ schematically diagrams interactions between three interacting GTPases and 
three phosphoinositides. I briefly outline the model equations describing these interactions. 
Further specifics can be found in [T3J [57] . Modifications of the model presented in those 
references, which are the subject of investigation here, are described in the main text. Each 
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GTPase undergoes conservative cycling between active membrane bound and inactive forms 
in the cell interior by (un)binding to the membrane. These dynamics are described by 

dG G'^ 

(34a) = Igt^ - SgG + D^G:,,, 

-g^ = ^^G^ + SgG + DcG%^, 

where G = R, p,C represents the membrane bound form and represents an inactive 
cytosolic form . Phosphoinositides interconvert between three states through the hydrolysis 
/ phosphorylation activity of PI5K, PI3K, PTEN etc. which are not explicitly modeled. 
The GTPase activation rate functions encoding the interactions in Figure 15.2b are defined 

by 

(34b) Ig = , , ,n , Ib =[ Ibi ' "^'^ * 



l + {p/a,T) ' " \^ l + /2(p/a3)"; ' 

/ - II 

l + {R/a,r- 

Phosphoinositidc kinetics are modeled by linear and mass action kinetics 

(34c) 
^P^ 

—l = Ip,- 5p,P, + ^21^2 - .fpi5K{R, C, p)P, + DpP,,,, 

ot 
BP 

= -^-21^2 + /p/5A-(i?, G, p)Pi - fpi3KiR, p)P2 + fpTEN{R, C, p)P3 + DpP^,,, 
-gf- = fpi3K{R, C, p)P2 - fpTENiR, C, p)P3 + DpP^,,, 

with feedback terms 

/o^j\ J. kpi3K l.,R\ r kpi^K L R 

(34d) ^^^^^-^ 1 + ^^^'^-^V^R, 



kpTEN p 

J PTEN — i H 

2 V Pt 

See Table [5] for a base parameter set for this model. 
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